Unifying variational metliods for simulating quantum many-body systems 
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We introduce a unified formulation of variational methods for simulating ground state properties of 
quantum many-body systems. The key feature is a novel variational method over quantum circuits via in- 
finitesimal unitary transformations, inspired by flow equation methods. Variational classes are represented 
as efficiently contractible unitary networks, including the matrix-product states of density matrix renor- 
malization, multiscale entanglement renormalization (MERA) states, weighted graph states, and quantum 
cellular automata. In particular, this provides a tool for varying over classes of states, such as MERA, 
for which so far no efficient way of variation has been known. The scheme is flexible when it comes to 
hybridizing methods or formulating new ones. We demonstrate the functioning by numerical implemen- 
tations of MERA, matrix-product states, and a new variational set on benchmarks. 
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Quantum many-body systems pose some of the most dif- 
ficult challenges in modern physics, and many examples 
remain inaccessible to analysis. Of the many methods that 
have been devised as attempts to meet these challenges, 
one of the most successful has been the density matrix 
renormalization group (DMRG)[ih. The DMRG was orig- 
inally conceived as a numerical technique for iteratively 
constructing the ground or low-energy states of a Hamil- 
tonian, so that the system's Hilbert space is truncated and 
the difficulties associated with exponentially increasing di- 
mension are avoided. A more recent interpretation of the 
DMRG has cast it as a variational method over matrix prod- 
uct states flSBEiS], and this shift in emphasis has stim- 
ulated much work on extending its applicability. 

Matrix product states are expected to provide good ap- 
proximations to the ground states of one-dimensional non- 
critical systems 10], however in other cases it is expected 
that alternative variational sets will be required. Motivated 
by this, new classes have been introduced, such as pro- 
jected entangled pair states [4] and weighted graph states 
list] for higher-dimensional lattices, while multi-scale en- 
tanglement renormalization (MERA) Jio!] and contractor 
renormalization 111 111 may be more appropriate for critical 
systems. At first sight it may appear that these numerical 
approaches to quantum many-body systems have little in 
common with each other. Moreover, the specification of 
a variational class is only a first step - we also require an 
effective method of finding the best description of the sys- 
tem's ground or low-energy states within that class. 

In this letter, we provide a unifying picture for several 
of these variational methods: We show that by (i) recast- 
ing state classes as quantum circuit classes (unitary net- 
works) one can formulate a (ii) general purpose variational 
method, related to the framework of flow equations [12]. 
We shall see that to provide a working numerical method, it 
is sufficient that the propagation of each Hamiltonian term 
can be efficiently tracked on a classical computer. The con- 
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FIG. 1 : (a) An example of what is here called a staircase circuit. 
The upper system here is d-dimensional. Further examples of 
unitary networks considered here include: (b) Unitary network 
of a weighted graph state up to local deformations, (c) Quan- 
tum cellular automata with Margolus partitioning, (d) Quantum 
circuit of MERA including disentangling operations, (e) A new 
variant combining (c) and (d) {"extended MERA"). 



tractibility properties of the state classes mentioned above 
translate immediately into an analogous property on their 
corresponding circuits. Our approach - a flow equation 
approach to variational simulations - may be regarded as 
an optimal control approach [ISl to varying efficiently 
contractible networks describing variational states of quan- 
tum many-body systems. It is flexible enough to hybridize 
known methods or to construct new ones, and provides a 
first efficient way of variation over MERA. 

Variational sets. - We begin with the casting of varia- 
tional sets as unitary networks, which provides the basis 
for the flow equations approach. Given spins (and pos- 
sibly ancillary systems) consider a family of states Sd = 
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{?7|0) : U G Ud}, where U ^ Udi^ set of unitary net- 
works characterized by some refinement parameter d, and 
|0) denotes the state vector with all spins down. The refine- 
ment parameter plays the role, e.g., of the auxiliary dimen- 
sion in matrix-product states (MPS). These networks con- 
sisting of gates U = rij^i ^^^^ to satisfy the condition 
that correlators of the form {0\WOi . . . OkU\0) can be ef- 
ficiently computed for any k, that is with effort polynomial 
in N and d. Before turning to the variational method we 
will discuss a number of important examples of states that 
can be discussed in this framework. We begin by consider- 
ing the matrix-product states of DMRG. For a system con- 
sisting of N spins, we introduce a d-dimensional ancilla 
system, and consider a circuit consisting of N gates Uj 
which act on both the ancilla system and qubit j, giving a 
''staircase" form, see Fig.[Ha). By projecting the output of 
such circuits onto some basis state of the ancillary system 
or by disentangling it, we obtain any matrix-product state 
on the spins. A second important example is the MERA 
class, see Fig.fHd). This circuit is arranged in a tree struc- 
ture with log N distinct layers, each of which introduces 
new spins into the circuit via two sets of gates known as 
isometrics and disentanglers ifToh . An analogous circuit is 
also possible using 2D binary trees. Further examples of 
quantum circuit classes are weighted graph states, where 
the refinement parameter d is defined by the non-zero en- 
tries of the adjacency matrix of the weighted graph, and 
quantum cellular automata [[l4ll . the finite depth d being 
the refinement parameter, and new variants as depicted in 
Fig.m 

Flow equations as a unifying method of variation. - 
Before we introduce the method of variation, let us first 
remind ourselves of flow equation ideas. Consider a 
continuous transformation of an initial Hamiltonian H 
H{t) = U{tyHU{t), where U{t) is defined via a Her- 
mitian generator as the time-ordered integral U{t) — 

T exp(— z G{s)ds). The derivative of H{t) is given by 
dtH — —i[G{t)^H{t)\. A familiar example from opti- 
mal control theory chooses G{t) — i[K^ H(t)], where K 
is a real diagonal matrix with unique entries. In this case 
H(t) will converge to a diagonal matrix of eigenvalues as 
t ^ oc, with the columns U{t) the corresponding eigen- 
vectors. This is often referred to as double -bracket flow 
1120. A straightforward application to quantum many-body 
systems is impractical, as the flow will in general transform 
the Hamiltonian into one having exponentially many terms. 
The key to these methods is to truncate the resulting sys- 
tems of differential equations in a perturbative fashion that 
is a good approximation for small perturbations. 

However, we are not aiming for approximate analyti- 
cal expressions here. Consider a quantum circuit as be 
a sequence of M gates Uj{t), each of which is continu- 
ously parameterized with infinitesimal generator Gj (t) be- 
ginning with some arbitrary Uj{0). Write the overall uni- 
tary implemented by the circuit as U{t) = 11^=1 



and consider the expectation of some many-body Hamil- 
tonian E{t) = {0\U{tyHU{t)\0). A circuit class is 
defined here by a specification of the locations of each 
gate Uj{t), and the best approximation to a ground state 
within a given class is the circuit that minimizes the ex- 
pectation E{t). Within the framework of flow equations, 
we will show how one can choose optimal generators 
Gj (t) for each gate. Differentiating the expectation we get 
dtE = 2Dlc{0\WHdtU\0), and our first goal is to min- 
imize this derivative subject to the Hilbert-Schmidt con- 
straints tr[Gj{tyGj{t)] = £, i.e., the generators remain 

U7^,Ujit),v,c find 



"infinitesimal". For U{t) 



M . M X . j X 

dtU = -iY,l n UkjcAHUkl (1) 



Note the convention whereby 11^=1 ordered as 

Um • • • U2U1, and the other way round for flj^M ^j- 
can substitute this back and minimize on a term-by-term 
basis at each point t of the flow. Let {B^} be an appropri- 
ate orthonormal Hermitian operator basis, and expand the 
j-th generator as Gj{t) = 9j,bB^^ with 5 real. Now 
define, for the given Hamiltonian H, 

T^^t,{t) = {0\U^H(ll Uk^B'fflUk^O). (2) 

Each term of the derivative with this parametrization is 
[dtE]j = 2 J2b 9j,b^^[^j,b{t)]^ and the constraints of the 
minimization problem are ^^^^^ = e. The Lagrange 
multiplier condition for a minimum is then simply 5 = 
— 29^e[rj,5(t)]/A. The following method of evaluating the 
optimal generator avoids calculating the quantities F^ 5(t) 
for each basis element B^: Writing Gj out in its opera- 
tor basis we have Gj{t) = -(2/A) J2b i^jA^)] 
Recall that the real part appears because we are taking an 
expectation of an operator and its Hermitian conjugate. It 
will be convenient to set Gj = -(2/A)(F^ + ^/ )' with 

3 M 

F,=ir[{j{U,)\0){0\U^H[ n Uu)b'^B], 



which can after some steps be written in terms of a partial 
trace over the subsystems Rj not acted on by the gate Uj 



3 M 

F3-^^R\{\[Uk)\0){0\U^H(^ f] Uu 



(3) 



The utility of this expression depends on whether or not 
we are able to efficiently calculate the partial trace, which 



depends on the structure of the circuit class Hj^i and 
on the Hamiltonian H. Accordingly, we now introduce 
contraction techniques able to cope with such expressions 
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FIG. 2: Contraction rules that are used in the procedure of finding 
the optimal generators (a-c), and a circuit correponding to Eq. ^ 
for MERA (d). The red and blue shading indicate the two causal 
cones that are encountered when evaluating the partial trace. 



Contraction. - We begin by reviewing some of the basic 
ideas of contraction and introduce a language for describ- 
ing them. What we here call a standard contraction prob- 
lem calls for the evaluation of an expectation of the form 
{0\WAU\0), U being a quantum circuit Hfii Uj actms 
on N spins (and any ancillae), and A some observable. A 
contraction procedure is a sequence of maps that construct 
operators Aq^ Ai^ . . . ^ Al with Aq = A and 

A = ((O|gj0i) ( n u)a-, (U^) (|o)qz®i). 



jeSi 



jeSi 



Here Qi and Si denote subsystems and subcircuits respec- 
tively, and are chosen so that at the final step we obtain 
{0\Al\0) = {0\WAU\0). The contraction is said to be 
efficient if the dimensions of the operators H/cg^ scale 
at most polynomially in the number of spins N. 

The key point here is that we may evaluate such an 
expectation (or similarly a trace or partial trace) without 
ever having to deal with the overall unitary U, whose di- 
mension is in general exponential in A^. As an example, 
suppose we have a two-body Hamiltonian term Hjj^i, 
and wish to evaluate {0\U^ Hk^k^iU\0) for U a staircase 
circuit. Then we set Aq — i?/e,/e+i and iterate Ai = 
(1 {0\k-i) Ul_iAi_^Uk_i (1 \0)k-i). so Qi here is 
simply the Ith spin, and Si contains only the gate Ui. The 
final operator A^ so obtained acts only on the ancilla sys- 
tem, and the desired expectation is thus (0|Ax,|0). At no 
stage in the procedure are we required to manipulate op- 
erators of dimension greater than 2d x 2d. To be entirely 
clear, a representative of each of these steps is depicted in 
Fig. [21 (a). A second example is provided by MERA cir- 
cuits [10], which require that we manipulate operators of 
dimension at most 64 x 64. Here the sets Qz, Si are de- 
fined with respect to levels of the MERA circuit and the 



causal cone of the given Hamiltonian term. A first such 
step is represented in Fig.[2](b). After constructing Ai, we 
can move on to the next MERA layer (Fig.[2](c)). 

We now describe the new contraction procedure used 
to evaluate the more general expressions of Eq. (O, see 
Fig. [21 There are two sets of gates highlighted in these 
circuits, which will be dealt with in two separate con- 
traction sequences, and we shall refer to these sets as the 
red cone (originating from the Hamiltonian as explained 
above) and blue cone (originating from the generator) re- 
spectively. The unhighlighted gates are simply cancelled as 
UlUk = 1. We can hence restrict ourselves to the causal 
cones. We shall also refer to the gate Uj whose generator 
we are calculating as the generate e. 

The first contraction sequence involves the gates in the 
red cone. Here we set Aq = and proceed via a 

sequence of partial expectations, with the same sets Qi, 
Si used in a standard contraction (see again Fig. [21 (b,c)). 
This continues until we reach the set Sq containing the 
generatee Uj . At this point we are unable to continue as 
some or all of the gates in Sq will have been cycled to the 
right-hand side (as we are calculating a partial trace). The 
final operator of this contraction Aq is then constructed by 
conjugating the previous Aq-i by those gates that have not 
been cycled. The second contraction focuses on the blue 
cone, and begins by initializing Bq = |0) (0|. acting on the 
subsystem as the generatee. The contraction then iterates 
in the reverse order to the standard contraction. 
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the primed Q'^ , S'j indicating the reversed order. This con- 
traction also continues until it reaches the set Sq, at which 
point Bq is constructed by (anti) conjugating Bq-i with 
the gates from Sq that have been cycled. The operator Fj 
is then given by Fj = tr^^ [BlAl], with Rj as for Eq. ([3]). 
For clarity, a MERA procedure is shown in Fig. [21(d). 

If the standard contraction procedure is efficient, then 
this modified procedure will also be efficient as the largest 
operators we must manipulate are defined by the same sets 
Si. For example, determining the optimal generator for a 
gate in a staircase circuit acting on N spins with ancilla di- 
mension d requires a time 0{Nd^), while for MERA the 
time required is 0{N log N). The above methods can be 
readily applied to 2D settings of MERA [10], where one 
has, e.g., layers of Margolus partitionings as in a quantum 
cellular automaton il4ll , with a tree-like reduction of the 
number of sites in every second step; then again, the con- 
traction of the two cones can be done efficiently. 

Implementation. - We now have the main ingredients for 
an actual algorithm. Fig. [51 illustrates example implemen- 
tations for the Heisenberg and critical Ising Hamiltonians 
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FIG. 3: (a) Staircase flow for the 30-qubit Heisenberg chain and 
(b) 28-qubit ring, in the number of flow steps. Plot (a) shows the 
decrease of expectation with the flow for d = 10 (dashed), 20 
(sohd), 40 (dotted). Exact values for the two lowest eigenvalues 
(for N = 30) are also shown, (c) MERA flow for the 32-qubit 
Ising chain, (d) MERA flow for the 32-qubit Ising ring and ex- 
tended MERA. 



as benchmarks, (a) and (b) illustrate an implementation 
for staircase circuits for the Heisenberg chain (with both 
open and closed boundary conditions) chosen for a first 
implementation as the corresponding ALPS-DMRG pro- 
vides good benchmark. Each step of the flow requires 
the calculation and application of the optimal generator for 
each gate. We see that for open boundary conditions the 
staircase achieves, for the energy E, the same accuracy of 
A = (E - Eo)/Eo as the benchmark ALPS-DMRG to 
six significant digits, and no problems with local optima 
have been observed [16] . (c) presents a MERA implemen- 
tation (representative when random initial conditions are 
drawn) for the 32-qubit critical Ising model for bond di- 
mension 2. Even for this small bond dimension, the rela- 
tive error of A = 4.4696 x 10~^ is achieved (note that this 
involves merely 61 unitaries acting on two spins, which 
is comparable in accuracy to DMRG for a dimension d 
defining MPS being described by an order of magnitude 
of more real parameters), for the ring A = 1.2901 x 10~^ . 
Similar performance is found for a 64 spin model, (d) For 
the extended MERA we find comparable performance but 
quicker convergence. We have also systematically com- 
pared the achievable accuracy for the Heisenberg model 
with MERA with bond dimension 2 (for which MERA per- 
forms slightly worse) with the one extended MERA where 
one appends an additional single layer of a quantum cellu- 
lar automaton: This leads in instances to a significant im- 
provement (of the order of (£"1 — Eq)/Eq in this model, 
critical in the thermodynamical limit), but, first and fore- 
most, shows the flexibility of the approach lisll . 

Conclusions and further work. - We have shown how 
ideas from flow equations may be adapted to provide varia- 
tional methods for approximating ground states of quantum 
many-body systems. The appeal of this approach is its flex- 
ibility as it is able to unify several existing ansatz classes 
within a single framework with a universal variational tech- 
nique. Recent work into the dissemination of correlations 
in quantum many-body systems has stimulated much work 
on the construction of suitable variational classes. It is 



hoped that the methods presented here will facilitate the 
systematic exploration of the potential of these approaches. 
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